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Empirical time series often contain observational noise. We investigate the effect of this noise 
on the estimated parameters of models fitted to the data. For data of physiological tremor, i.e. a 
small amplitude oscillation of the outstretched hand of healthy subjects, we compare the results 
for a linear model that explicitly includes additional observational noise to one that ignores this 
noise. We discuss problems and possible solutions for nonlinear deterministic as well as nonlinear 
stochastic processes. Especially we discuss the state space model applicable for modeling noisy 
stochastic systems and Bock's algorithm capable for modeling noisy deterministic systems. 



I. INTRODUCTION 



One of the aims in time series analysis of complex systems is to find a dynamical model for the data |l9| . A well 
fitting model might yield insight into the underlying process. But also if prediction is the main goal, e.g. in order to 
discriminate chaos from stochasticity |6C],[l7]], or if fitted models are used to determine dynamical invariants of the 
system 37 Sfl, an optimal estimation of the parameters is desirable. Our own interest was stimulated by modeling 



physiological time series in order to gain insight into the underlying systems. 

The dynamical models can be (time-continuous) differential or (time-discrete) difference equations. Methods to 
estimate parameters in differential equations can be subdivided according to whether they require an estimation 
of the derivatives from the data (ll]]3(J] or not |7U23]|. Difference equation allow to use a great variety of different 
approaches to model the mapping from past values to the present one. These range from parametric ones |29|:pf via 
neural nets (ijj7^,^] , radial basis functions and nonparametric models |^38| to local linear models [|l7||. 

For all these methods a significant amount of observational noise can be a severe problem. Especially for difference 
equations the functional relation between past and present values will be underestimated if observational noise is 
not included in the model. In the first part of this paper we exemplify this for data of physiological human hand 
tremor. These data contain up to 50 % observational noise. Inspired by the physics of the process, linear stochastic 
autoregressive (AR) processes have already been suggested in 1973 to model these data |j50| . In [^J these data were 
analyzed by a linear state space model which takes the observational noise into account. For one data set we compare 
these two approaches in detail. In the second part of this paper we discuss problems introduced by observational 
noise and possible solutions for non-linear deterministic as well as non-linear stochastic processes. 



II. MODELING PHYSIOLOGICAL TREMOR 

The outstretched hand of a healthy subject exhibits a small amplitude oscillation called physiologic tremor |2^] . The 
movement of the hand can be measured by piezoresistive accelerometers attached at the hand. Simultaneously the 
flexor and extensor muscle activity is recorded. The spectra of the muscle activity data are often flat in physiological 
tremor reflecting an uncorrelated activity of motoneurons f34| . For an analysis of physiological tremor data in which 
the muscle activity is correlated and included in the modeling, see flog] . 

Fig. [I] displays a 2 s section of the whole 35 s measurement which is sampled with 300 Hz. The time series consist 
of 10240 data points covering approximately 250 periods of the oscillation. The data are normalized to zero mean 
and unit variance. Fig. || displays the periodogram, i.e. the absolute value of the Fourier transform squared. For 
the treatment of special problems of estimating the spectrum from the periodogram for tremor time series, see JS3| . 
The periodogram shows on the one hand the high amount of observational noise in these data. Roughly estimating 
the variance of the signal by summing up the periodogram in the range from 2.5 to 12.5 Hz, and that of the noise 
from the remaining frequencies results in a signal-to-noise ratio of 0.93 in relative amplitudes. On the other hand 
the periodogram shows one broad peak around a frequency of 7.5 Hz. This peak is usually explained as a resonance 
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phenomenon |p8[ . The outstretched hand is a damped oscillator which is excited by the uncorrelated muscle activity. 
Usually no higher harmonics show up in the periodogram indicating that the process is linear. 




FIG. 1. A 2 s section of human physiological hand tremor data. 
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FIG. 2. Periodogram of the physiological tremor data set. 



The autoregressive ( AR) processes invented by Yule 1 72 are expected to fit such data well. An AR process of order 
p is given by : 



= Y^aix{t-i) + e(t) , 



(1) 



where e(t) denotes a uncorrelated Gaussian distributed random variable with mean zero and variance a 2 . Such 
a process can be interpreted, depending on the chosen parameters, as a combination of relaxators and damped 
oscillators Q. For example, an AR process of order 2 that corresponds to a damped oscillator with characteristic 
period T and relaxation time r given by: 
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The spectrum is given by : 
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1 - ax 
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(4) 



AR processes can be generalized to the autoregressive moving average (ARMA[p,q]) processes by including past q 
driving noise terms in the dynamics. For theoretical reasons ARMA[p,p-l] should be preferred to AR[p] processes for 
modeling of sampled continuous-time processes 46 . Experience shows that differences in the results are small. 



A more substantial generalization is the linear state space model (LSSM) 34 which enables the explicit modeling 
of observational noise rj{t) that contributes to the measured y(t) : 



x{t)=Ax{t-l) + e{t), e(t) e Af{0, Q) 
y(t) = Cx*{t)+r)(t), V [t) e JV(0, R) 



(5) 
(6) 



Eq. (JsJ) describes the linear dynamics. Eq. (|^) maps the dynamics to the observation including the observational noise 
r)(t). Another advantage of the LSSM is its capability to model superpositions of linear processes. This is not possible 
for AR or ARMA processes. The spectrum of a LSSM is given by : 



S(lu) = C(l - Ae^j^Q ((1 - Ae iui y 
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(7) 



The superscript T denotes transposition. Spectra of AR or ARMA processes are special cases of Eq. (^). 

While parameter estimation in AR models e.g. by the Burg or Durbin-Levinson algorithm is well established, 
parameter estimation in the LSSM is more cumbersome, see |34| for a detailed description. Usually the Expectation 
- Maximization (EM) algorithm is applied f2"l|]54| . The EM - algorithm is a general iterative procedure to estimate 
parameters for models in which not all variables are observable, here x(t). Denoting the joint density of x(t) and y(t) 
given the parameters by p(x,y\&) and the density of x(t) given the data y(t) and the parameters by p(x\y, 0) 
the quantity : 



2.(6,0*) =< ln(p(x,y\Q)) >„( x | w ,e«) 



(8) 



is calculated in the ith Expectation step. In the Maximization step, £(0, l ) is maximized with respect to yielding 
I+1 . 

In the case of the LSSM, this means that starting from some initial values for the parameters A, Q, C, R the hidden 
dynamic variable x(t) is estimated by the Kalman filter J30| ] in the expectation step. Denoting the estimator of a 
quantity z(t) based on the data y(l), ...,y{t') by z t \ t >, the covariance matrix of the estimated x(t) by O t u/ and the 
variance of the prediction errors (y(t) — yt\t') by A t i t / the Kalman Filter reads: 
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Since the model is gaussian and linear, the density p{x\y, 0) is completely specified by x t \ t and their covariances 
Qt\t- Note, that the covariances Q t u only depend on the model parameters, not on the data. An improvement of 
the estimated quantities by the so-called smoothing filter and the lengthly equations for the parameter update of 
A, Q, C, R in the Maximization step are given in [ p4[ . This procedure is iterated until some convergence criterion is 
fulfilled. 

There are different criteria to judge the goodness of fit of a dynamical model. 
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• Whiteness of the prediction errors. The model should explain all correlations in the data yielding uncorrelated 
prediction errors. By the Kolmogorov - Smirnov test ]48| it can be tested whether the periodogram of the 
prediction residuals is consistent with white noise ]T^ |. 

• A knee in the variance of the prediction errors for ascending model order indicates the correct order. Criteria 
like AIC ||, BIC [Q, MDL |52| can be used to take into account the different number of parameters in the 
compared models. The application of these criteria is not without problems 0. If the model class is not correct, 
these criteria will choose "some" order. 

• The distribution of some feature can be derived from numerous realizations of the model and the compatibility 
of the value of the feature calculated from the data with this distribution can be examined, see e.g. |7l] ]. 

In the case of linear modeling there are two more criteria. 

• Since the parameters in linear models are related to relaxation times of the corresponding oscillators and 
relaxators, negligible relaxation times indicate a to large model order. 

• According to Eq. (0) the spectra of linear processes can be calculated from the parameters of the estimated 
models. Since the spectrum is the expectation of the periodogram, a comparison of the spectrum to the 
periodogram of the data serves as a qualitative criterion. 

We now present the results of the analysis of the data shown in Fig. [I]. We do not give the estimated parameter 
here since they are rather uninformative. Instead, Tab. 1 displays the resulting periods and relaxation times. These 
characteristic times can be calculated for any model order by formulas analogous to Eqs. (||,||) [ p4[ . 
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TAB. 1. Relaxation times and periods of the fitted AR and LSS models of ascending order for the physiological 
tremor data. "prob. WN" denotes the probability of the residuals to be consistent with white noise. 

For an estimated oscillator a relaxation time smaller than the period of the oscillation indicates that this oscillator 
is insignificant. Pure relaxators with an relaxation time of a few time steps have also to be regarded as negligible. 
For the LSSM, a model order of two with a period of approximately 40 time steps, corresponding to a frequency of 
7.5 Hz, and a relaxation time of 91 time steps is clearly identified. For the AR models no significant structure is 
detected. Tab. 1 also gives the probability of the residuals to be consistent with white noise. For the LSSM, models 
with an order larger than one yield prediction residuals which are consistent with white noise. The residuals of the 
AR processes are not consistent with white noise for all model orders investigated. 

Fig. U displays the variance of the prediction error in dependence on the model order. The graph for the LSSM 
exhibits a clear cut knee at an order of 2 whereas the graph of the AR model does not allow for a clear cut decision. 
The graph saturates around order 5 and decreases slowly for larger orders. We fitted AR model up to order 50 and 
applied the AIC. Even for this unrealistic high orders the function AIC(p) still decreases. 

Fig. |] shows a quantile-quantile - plot with respect to a Gaussian distribution of the normalized prediction residuals 
for the LSSM of order two and the expected straight line. Only ten of the 10240 points deviate from the straight line 
indicating Gaussianty of the prediction residuals. 

All these results indicate that the LSSM of order two describes the data adequately. Fig. [| displays the low frequency 
part of spectrum estimated from the parameters of the fitted LSSM of order two and the periodogram of the data. 
This confirms the above results. 
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FIG. 3. Residual variance of the fitted AR (+) and LSS (O) models in dependence on the model order. 
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FIG. 4. Quantile-quantile plot of the normalized prediction residuals of the second order LSSM with respect to a Gaussian 
distribution. The straight gives the expected behavior in the case of Gaussianity. 



In terms of physics and physiology these result show that the physiological tremor is a linear damped oscillation 
driven by uncorrelated muscle activity. The parameters of the model describe the mechanical properties of the 
muscle-hand system, i.e. its mass, stiffness and damping. 



III. MODELING NOISY TIME SERIES 



The behavior of the AR model in the previous section, i.e. not recovering the correct model order, is caused by the 
observational noise. Whereas the LSSM includes the observational noise explicitly in the model, the AR model assumes 
the data to be free from observational noise. We use a simple first order process to demonstrate the consequences. In 
an AR[1] model 

x{t) = ax(t - 1) + e(t) (16) 
the parameter a can be estimated without bias by : 

. ^x(t-l)x(t) 



5>(t - iMi - 1) 



(17) 
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FIG. 5. Periodogram of the data (dotted) and spectrum of the fitted second order linear state space model (solid). 

If the dynamics is covered by observational noise: 

y(t)=x(t)+r 1 (t), r){t)~Af(Q,R) , 
the expected value of a estimated analogously to Eq. ( p"7j ) from y(t) is 

< y(t - l)y(t) > 1 



< a >- 



< y(t - l)y(t - 1) > 
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Thus, the parameter a is underestimated. The degree of the underestimation depends on the signal-to-noise ratio. 
This effect is known from linear regression theory where it is called the " Error- in- variables" - problem |26|| and was 
first mentioned in the context of time series in |4^] . For the data analyzed in the previous section the variance of the 
signal nearly equals the variance of the noise. Thus, an underestimation of a by a factor of two is expected. Indeed, 
even for the misspecified first order models, the parameter estimated from the LSSM is 0.882 and that for the AR[1] 
model is 0.443. Note, that in the case of second order models the misestimation of the parameters for the AR[2] 
models leads to a detection of two relaxators, whereas the LSSM recognizes the oscillator, see Tab. 1. Thus the AR 
model leads to a wrong interpretation. 

The underestimation of the functional relation between past and present values in the case that observational 
noise is not taken into account carries over to more general models. Assuming stationarity of the underlying process 
systems can be linear or nonlinear, resp. deterministic or stochastic. Our approach to discriminate deterministic from 
stochastic systems, i.e. the presence of dynamical noise, is rather practical. For example, in the data analyzed in the 
previous section it is imaginable to model the muscle activity by a deterministic system, including a model for the 
brain and the spinal tract. This would result in a very high dimensional dynamical system whose output is equally 
well captured by white noise. Thus, whenever a huge amount of influences enters the system significantly a stochastic 
description might be recommended. But also if there exists no shadowing trajectory for a realization of deterministic 
system, the use of stochastic models might be indicated § • 

The dynamics can be described by differential equation or by maps. In the following we discuss possible treatments 
of observational noise in the different cases. If present the dynamical noise is assumed to be Gaussian and additive. 
By the latter assumption we exclude models which are able to describe fluctuations in the parameters, e.g. bilinear 
models 159]. 



A. Linear deterministic case 



The linear deterministic case is easy to treat since the solutions of the dynamical equations are explicitly known 
and can be fitted to the noisy data. Parameters of continuous- and discrete-time modeling are connected by equations 
like Eq. (H,||) . An application to physiological data is given in |32| . 
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B. Linear stochastic case 



The linear stochastic case is solved by the state space model discussed above. Because of the linearity and gaussianity 
of the model there is a direct connection between the parameters from discrete-time and those from the continuous- 
time models. Parameter estimation for the continuous-time case is presented in JSTj. An application of the state space 
model to astrophysical data is given in [|d| . 



C. Nonlinear deterministic case 



The effect of observational noise is different for fitting differential and difference equations to the data. 

For modeling data by differential equations two approaches can be distinguished depending on whether time deriva- 
tives of the process has be estimated from the data or not. Since estimating derivatives form the data amplifies the 
noise, the former method as applied in |18 l^^,^,|l|,^5) is vulnerable to significant amounts of observational noise 



as demonstrated in 

To fit differential equations to noisy deterministic data without estimating derivatives from the data, at least two 
different approaches exist. In the initial value approach one chooses some initial parameters in the model and an 
initial value for the trajectory, integrates the differential equation and tries to minimize the prediction error by some 
minimization algorithm p3| . Without many precautions, this procedure has been shown to be unstable, i.e. yielding 
divergent trajectories, and is susceptible to stopping in a local minimum [[5l|]66| 1. A more sophisticated algorithm was 
developed by Bock His elegant multiple shooting approach starts with a discontinuous trajectory which stays 

close to the data. The continuity of the underlying trajectory enters into the algorithm by a constraint in the cost 
function. This constraint is nonlinear in the parameters but enters the optimizing strategy only in a linearized way. 
Therefore, the trajectory is allowed to be discontinuous at the beginning of the iteration but is forced to become 
continuous in the end. 

We illustrate the behavior of Bock's algorithm by the restricted Lotka-Volterra system which has been used to 



compare different fitting procedures P,p9|,p7 23 



±i = k\X\ — k 2 x\x 2 (20) 
±2 = k 2 xix 2 - k 3 x 2 (21) 

with the parameters fci = 1.0, k 2 = 1.5 and £3 = 2.0 and initial values xi(0) = 0.3, £2(0) = 0.8. The system is 
integrated by the routine BSSTEP from for 20 s and sampled with At = 0.2 s. The standard deviation of the 
data is 1.09. Observational noise with standard deviation of 0.5 was added to the signal. To fit the parameters, 
only the first component was used. All starting points for the second component were set to 0.5. Initial values 
for the parameter estimates were k\ = 0.5, k 2 = 0.75 and k 3 = 1.0. 30 starting points were used for the initially 
discontinuous trial trajectory. Fig. |^ shows initial situation (A), the third iteration (B) and the converged solution 
(C) after 16 iterations of the algorithm. The true trajectory is reproduced well. With respect to the confidence 
regions, the estimated parameters are compatible with the true parameters: k\ = 1.01 ± 0.14, k 2 = 1.53 ± 0.34 and 
k 3 = 2.02 ±0.45. 

Simulation studies comparing both algorithms show that Bock's algorithm is more stable and converges to the 
global minimum for a larger set of initial guesses for the parameters than the initial value approach ]5l]j6q| . Note, 
that for both approaches it is not necessary to reconstruct the dynamics by delay-embedding [ |S~3~| ] avoiding all problems 
related to choosing the delay time (43| . 

In a simulation study Bock's algorithm has been successfully applied to model noisy chaotic data from dissipative 
as well as conservative systems ||,[38) . Applications of this method to chemical and physiological data are reported 
in |]J|. The basic idea of Bock's algorithm, i.e. taking advantage of the fact that the process under investigation 
has produced a continuous trajectory, carries over to modeling noisy deterministic systems by maps. Nevertheless 
this idea has not attracted much attention and modeling data by maps is generally treated as a regression not as a 
dynamical problem. 

Stimulated by JTs| ] a large number of different approaches have been suggested to model nonlinear dynamics by maps. 
They are distinguished by the type of basis function used to approximate the nonlinear functional relationship between 
past values and the present one. Amongst others, polynomials p9| , sigmoidals |44|j7C| , ^9| , radial basis functions ^7]] 
and local linear models jlTj have been suggested. Analogous to the AR model discussed above all these methods arc 
subject to the " Error- in- variables" - problem. For the tremor data and the AR model of order two in Sec. |l| above 
this problem led to a detection of two relaxators instead of an oscillator. In [R7| it was reported that observational 
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noise covering chaotic data led to a detection of a periodic process. That might be caused by the same phenomenon 
of underestimated parameters. 

The generalization of the theory of " Error- in- variables" from nonlinear regression, see ]l5j for a recent review, to 
modeling of noisy nonlinear deterministic time series by maps and an application to empirical data are presented in 
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FIG. 6. Convergence of Bock's algorithm for the noisy Lotka-Volterra system, 
iteration. C: After convergence, dashed line: true trajectory. 
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D. Nonlinear stochastic case 



For the nonlinear stochastic case, including a possible nonlinear mapping g(.) of the state vector x(t) to the 
observation y(t) : 



m = /(*(*)) 

y(t) = g(x(t)) 



e(t) e Af(0, Q) 



(22) 
(23) 



a generalization of the Kalman filter Eq. (ffllq) appears to be the solution. Therefore, for the discrete-time version, 
the matrices A and C are replaced by linear approximations of f(x(i)) and g(x(t)) p8 45|: 



A 

C 



dx 
dg_ 

dx 



X = X t _ 1 l t _ 1 



(24) 
(25) 



This solution is, in general, not valid for two reasons. First, as mentioned in Sec. ||, the variance f2 t | t of the predicted 
state variable x t \t depends only on the model, i.e. not on the amount of data, and might be so large that the 
linearization is not valid. Thus, the distribution of x t \t becomes nongaussian and a first and second order description 
is not sufficient. Second, the noise will become nongaussian in time-discrete nonlinear models. This results from the 
theory of integrating stochastic differential equations, see [ flO| for a recent review. One main result of this theory is 
that unlike deterministic differential equations, stochastic ones can not be integrated by arbitrary high order methods, 
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e.g. by integrating the deterministic part by a Runge-Kutta method of 4.th order and adding some noise to the resulting 
value. Due to complicated stochastic integrals in the Taylor expansion for higher order integration schemes, only first 
or second order methods are feasible. Therefore, usually, the process has to be integrated on a much smaller time 
scale than it is observed. This turns the effective noise on the time scale of the observation nongaussian. Because of 
this effect methods like those proposed in PJiofl are applicable in special cases only. 

The nongaussianity prevents a maximum likelihood estimation of the parameter. A least square approach yields 
biased estimates 1 33 . In general, for nongaussian distributions all higher moments have to be taken into account. 
A truncated expansion in moments was suggested in p8|| , the application of the Gibbs sampler using a mixture of 
normals to approximate the distributions was introduced in |14j . To our knowledge no application of these methods 
to empirical data have been reported. 



All the map-based methods mentioned in Sec. Ill C can be applied to nonlinear stochastic system by allowing for a 



residual error in the prediction. Here, the above mentioned non-gaussianity is also a problem since commonly applied 
least square fitting is no longer maximum likelihood. 

An application of the "Error-in- variables" -theory to this double stochastic situation is not straightforward. To 
apply this theory, the ratio of the variance of the noise of past values to that of the present value must be known. The 
noise covering the past values is simply the observational noise, but the noise on the present value is, in contrast to 
pure regression, the sum of the observational and the dynamical noise. Since spectra of most (observational noise free) 
processes decay for high frequencies |56| , |55| ], the variance of the observational noise can be estimated from this region. 
But the variance of the dynamical noise is usually obtained after fitting the model since minimizing this variance is 
the optimizing criterion. 

The severity of the effects of the problems discussed depends on the system, i.e. the type of the nonlinearity and 
the variances of dynamical and observational noise. We demonstrate this by briefly reporting the results of modeling 
a time series of Parkinsonian tremor. This type of tremor is believed to be a nonlinear process [^7[62|. The time 
series was recorded under the same conditions as the physiological tremor data analyzed in Sec. 0. Fortunately, the 
observational noise in these data is negligible small. Using a large variety of ansatzes for the right hand side of the 
differential equation, we did not succeed in fitting a nonlinear deterministic differential equation by Bock's algorithm 
to the data. Therefore, we assumed that the process is stochastic and tried global polynomials that are orthogonal 
with respect to the measure given by the data. These polynomials were discussed in the context of regression theory 
by Forsythe p5j . He also introduced a recursive algorithm to determine the polynomials and their coefficients. In 
p9| , ^3| j2j| they were made popular for the analysis of time series. As discussed in || the significance of each single 
polynomial can be judged independently from the other polynomials. 

Using models of polynomial order 4 and an embedding delay r of 10 time steps we fitted models of ascending 
autoregression order to the data. We found a model of regression order 3 containing 16 significant parameters which 
describes the data adequately. The goodness of fit was judged by a clear cut knee in the residual variance (Fig. |t]), 
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FIG. 7. Variance of the prediction errors for 4.th order polynomial models of ascending autoregression order, 
the whiteness of the residuals (Fig. ||) and a comparison of the spectra of the empirical data and data realized from 
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FIG. 8. Kolmogorov-Smirnov test statistic D for the consistency of the prediction errors with white noise of 4.th order 
polynomials models of ascending autoregression order. The 5 % significance level is marked. 

the model (Fig. |9|). Note, that in Fig. || also the small peak near 13.5 Hz which is only due to aliasing because of 
the effective downsampling by choosing t = 10, is reproduced by the model. 
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FIG. 9. Spectrum estimated from the original data (solid line) and from a realization of the fitted model (dotted line). 



Unfortunately we did not succeed in assigning a physiologic meaning to the fitted parameters. 



IV. DISCUSSION 



The project "Equations of Motion from a Data Series" [|19| has attracted much attention. Our own interest was 
to obtain models for physiological data like EEG or tremor. The hope was to learn something about the underlying 
systems from the fitted models. In the first part of this paper we gave an example where this idea could be fulfilled 
and the parameter of the fitted model could be interpreted in terms of physics and physiology. This was facilitated 
by the fact that the well understood theory of linear systems was applicable. 



10 



In the general case of nonlinear deterministic or even nonlinear stochastic systems it appears to be much harder 
to obtain such a result. This is on the one hand caused by the fact that the treatment of observational noise is not 
solved in general. On the other hand without much knowledge about the system under investigation, it is hard to 
decide on a certain interpretable ansatz for the parameterization of the nonlinear dynamics. 

Often an interpretable model is not the goal of modeling, e.g. if only prediction as for stock market data or the 
temporal evolution of the prediction error as for discriminating deterministic from stochastic behavior is the aim. But 
here also observational noise is a problem, since the functional relationship will be underestimated if it is modeled by 
maps. A model taking this noise into account would lead to a better model and a smaller prediction error. 

Results for fitted differential equations might be compared to known systems in order to understand the underlying 
mechanisms. For maps this seems to be much more difficult, especially when sigmoidals, radial basis functions or 
local linear models are used as basis functions. But even if global polynomials are used an assignment of a physical 
meaning to the parameters is difficult since a single parameter in the underlying differential equation shows up in the 
coefficients of more than one basis function. Furthermore, the model structure depends on the sampling time. 

For many methods applied to model time series, observational noise causes problems since these methods treat the 
data as in regression. For deterministic systems Bock's algorithm provides an elegant alternative which explores the 
additional information that the process under investigation produced a continuous trajectory. 

There is growing evidence that many physiological processes are neither linear stochastic processes since surrogate 
data testing reject this hypothesis nor nonlinear deterministic processes since low dimensional attractors can not be 
established. Therefore, we believe that methods to model (noisy) nonlinear stochastic processes are worth being 
studied in more detail. 



DATA AVAILABILITY 

The tremor data and the simulated data of the Lotka- Volterra system are accessible at: 
http://phyml.physik.uni-freiburg.de/^jeti/tremordata/ 
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